Particles, holes and solitons: a matrix product state approach 
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We introduce a variational method for calculating dispersion relations of translation invariant 
(l+l)-dimensional quantum field theories. The method is based on continuous matrix product 
states and can be implemented efficiently. We study the Lieb-Liniger model as a benchmark where, 
despite criticality, excellent agreement with the exact solution is found, including, clear solitonic 
effects in Lieb's Type II excitation. In addition, a non-integrable model is studied where a (7(1)- 
symmetry breaking term is added to the Lieb-Liniger Hamiltonian. For this model we find evidence 
of a non-trivial bound-state excitation in the dispersion relation. 
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The last decades have witnessed an explosion in the ex- 
perimental realization of strongly correlated one-dimensional 
quantum systems [l|]. Often a discretised description in terms 
of a lattice Hamiltonian is possible, which can then be stud- 
ied using White's density matrix renormalization group 
The underlying variational ansatz, the set of matrix product 
states (MPS) [3t3], explains both the success of this method 
and has been used to develop generalizations beyond the set- 
ting of ground states, e.g . to the study of time evolution@j9| , 
dissipative dynamics [ly, [ll[ and dispersion relations [l4.lla|. 

The use of elongated optical or magnetic atom traps have 
opened the possibility of creating one-dimensional quantum 
gasses in the lab [lj-[l3| . It is natural to study these systems 
directly using quantum fields, without resorting to a lattice 
discretisation. A continuum limit of matrix product states, 
known as continuous matrix product states (cMPS) [20l - [2^ . 
was recently developed and has demonstrated to be an effi- 
cient description of the ground-state properties of the Lieb- 
Liniger (LL) model [S^l which describes a one-dimensional 
interacting Bose gas. 

Apart from ground-state properties, there has also been 
experimental interest in localized excitations in these systems 
[2J-[23|. While Lieb determined the spectrum of excitations 
for the LL model [2^, a systematic method for studying ex- 
citations of non-integrable quantum fields is still lacking. In 
this Letter we fill the gap by extending the recently introduced 
ansatz for excitations of translation invariant spin chains in 
the thermodynamic limit [l3| to the setting of cMPS. This 
yields a new variational ansatz for elementary excitations of 
translation invariant quantum fields and allows us to simulate 
dispersion relations. The corresponding variational states are 
faithful eigenstates and have therefore in principle an infinite 
lifetime, out of which we can of course construct localized 
wavepackets by taking linear combinations. In particular, we 
are able to reproduce the spectrum of the LL Hamiltonian and 
also study the effect of adding non-integrable perturbations. 

A cMPS for a translation invariant infinite system with 
open boundary conditions is defined as [23l 



\^{Q,R))=vl (re 
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where Q,R e C^*"^, vl and vr are D-dimensional bound- 
ary vectors acting on an ancillary system, |f2) is the Fock 
vacuum and V the path-ordering operator. For bosonic sys- 
tems the field operators satisfy the commutation relation 
[ip{x),-ip^ (y)] = 5{x — y). For a generic normalizable cMPS, all 



eigenvalues of the transfer matrix T = Q ®1 -'r'iL® Q + R® R 
have non-positive real part and there is a non-degenerate zero 
eigenvalue. The corresponding left and right zero eigenvec- 
tors, {l\T = 0, T\r) = may be reshaped to give posi- 
tive, hermitian matrices I and r which have full rank and 
are normalized so that {l\r) — Tr{lr) — 1. Since the bound- 
ary vectors have no variational importance they are chosen 
so that the state has norm 1, i.e. {'i'{Q, R)\'i/{Q, R)) = 
{vl(Svl)\r){l\ {vR(g)VR) = {vlrvL){vlilvB^) = 1. 

Suppose we have found the ground state (or its best approx- 
imation) of a system as a cMPS parametrized by the matrices 
(Q,R). An ansatz for a particle-like eigenstate or excitation 
is obtained by locally replacing the matrices Q and _R by V 
and W — which has the effect of perturbing the ground state 
in a spatial region of the size of the correlation length — and 
then building a state with definite momentum via a plane 
wave superposition: 
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X {y ®t + W ®i}\x)^jj2{x,oo)vR\Q.) , 

where Ua{x,y) = P exp(/J dz (Qa ® t + Ra ® ■i\z))) for 
a — 1,2. The crucial feature here is a single localized distur- 
bance described by the matrices V,W £ C^^^, characteriz- 
ing it as a single-particle ansatz. The states |$p(y, VK)) are 
momentum eigenstates and obey a (5-orthogonality. Note that 
they depend linearly on the variational parameters V and W, 
so they span a linear subspace of Hilbert space. Asymptoti- 
cally the states \$p{V,W)) look like |*(Qi,i?i)) at a; = -oo 
and \^{Q2, R2)) at a: = +00, which are supposed to be equally 
good but potentially different ground states (in the case of 
symmetry breaking). Our ansatz thus includes the possibility 
of capturing topologically non-trivial excitations, henceforth 
referred to as topological excitations for the sake of brevity. 
For |$(Qi, Ri)) = \^(Q2, R2)), the ansatz describes (topolog- 
ically) trivial excitations and can be interpreted as a state in 
the momentum sector p of the tangent space obtained by an 
infinitesimal position dependent variation of |^(Q,-R)) [22 ]. 

To compute excitations we apply the Rayleigh-Ritz 
method, which minimizes the energy functional of a Hamilto- 
nian H within a given linear variational subspace and results 
in a generalized eigenvalue problem for an effective Hamilto- 
nian. In our case the generalized eigenvalue equation is given 
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by 



with 



= EN^ 



[ ^- ] being the 2D -dimensional vector corresponding to 
V and W , E is the energy, Hp the efTective Hamiltonian and 

l2 r, n2 



dimensional matrices defined by 

{%{V,W)\H ~ Eo\%,{V,W)) = 2tt5{p^p') [ytvyt] H, 
{^j,iV,W)\%,iV,W)) = 2-k5{p ~ p') [yt^/t] 



where Eq is the ground state energy obtained with the ground 
state approximation |5'(Q,_R)). 

Note that Hp and A'p have zero eigenvalues correspond- 
ing to a redundancy in the representation of the states 
|<I>p(Vi W)), wliicli can be traced back to the gauge invariance 
of the states \'^(Q,R)) under the transformation Q g~^Qg 
and R ^ g'^Rg [3- One finds that for all X £ C^^-°, 
\%i[X,Q] + ipX,[X,R])) = 0. Hence, \%{V',W')) = 
\$p{V,W)) if V" = V + [X,Q]+ipX and W = W+[X,R]. If, 
for p = 0, one also restricts to states orthogonal to the ground 
state by imposing {l\ (V ® t + W (E) R) \r} = 0, there are 
gauge degrees of freedom for every momentum p. These can 
be fixed by constraining V and W to satisfy 



{1\{V (g)l + W (g)R) = ^ 



(1) 



with {l\ the left eigenvector of Tn, the transfer matrix corre- 
sponding to |*I'(Qi, _Ri)). We can parameterize V and W in 
terms of a single matrix F e C^''^ as = r^/^Fr"^/^ and 
V = —1~^RIiW, where r corresponds to the right eigenvector 
of T22, the transfer matrix of \'ii{Q2, R2)) ■ It can be shown 
that for this parameterization the effective norm matrix A'p is 
the identity [2^ and the Rayleigh-Ritz method is transformed 
into an ordinary eigenvalue problem. 

Evaluating the effective Hamiltonian Hp is more involved 
and requires the specification of a Hamiltonian H . We show 
the explicit result for the repulsive LL Hamiltonian [S^] 
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where c > is the interaction strength, v < the potential 
and where we have omitted the argument (x) of the field op- 
erators for the sake of brevity. The Hamiltonian is gapless 
and only depends on a single parameter 7 = c/p with p the 
ground state density. To find the ground state of the system 
we apply the time dependent variational principle (TDVP) 
to cMPS, yielding an optimal Q and R. To calculate the ex- 
pectation value of Hll with respect to the excitation ansatz 
|$p(Vi W)) it is useful to note that 

ip{x)\^p{V,W)} = vU f dye''"'Ui{-oo,y){V + W (Sip^ 

J —00 

xU2(y,x)R2U2ix,oo)+ / dye''"'Ui(-oo,x)RiUiix,y) 

J X 

X {V ®t + W ®'ip'^)U2{y,oo) 

+ e''"'Ui{-oo,x)WU2{x, oo)\^VR\n). 

This result can then be used to compute '0(i/)-i/)(x)|$p(y, W)) 
and |^(a:)|cE>p(V, W)). The expectation values of all three 



operators in can be found by calculating the norm of 

these three states integrated over x and result in the following 
expression, where the constraint in Eq. [T]is assumed to hold; 

($p(F,W)|i?LL -£o|$p'(V,H/)) = 2'k5{p' -p) X 

(/| ([Qi, -Ri] ® [<3i,-Ri] + cRI (»rI+ vRi (g) i-T)'' 

X l^{V 01 + W ® Ri) (T21 + ipf (1 ® F -f i?2 ® W) 
+ {1 (g)V + Ri ®W) (T12 - ipf {V(S)t + W(g)R2) 
+ W®Wj 

+ W ®W{-T)'' {[Q2, R2] (S IQ2, R2] + cRI ®RI + vR2 ® R2) 
+ [lQi,Ri](S [{QiW - WQ2) + {VR2 - RiV)] 

-ip[Qi,Ri](g)W + cRi(g) {WR2 + RiW) + vRi ® Wj 
X {T12 - ipf {V ®1 + W(» R2) 
+ { [{QiW - WQ2) + {VR2 - RiV)] ® [Qi,Ri] 

+ ipW® [Qi,Ri] +c{WR2 + RiW)®Ri + vW®Ri'^ 
■K {T21 + ipf {t ®V + R2 ®W) 
+ vW®W + c {WR2 + RiW)® (TFi?2 + RiW) 
+ ipW ® [{QiW - WQ2) + (Fi?2 - RiV)] 
- ip [{QiW - WQ2) + {VR2 - RiV)] ®W 
+ p''W®W 

+ [{QiW - WQ2) + {VR2 - RiV)] ® {QiW - WQ2) 



+ [{QiW - WQ2) + {VR2 - RiV)\ ® {VR2 - RiV) 



\r). 



This expression involves a generalized transfer operator Tab ~ 
Qa ® 1 + t + Ra ®Rb tor a,b = 1, 2. Additionally, we 

have introduced the notation T'' to denote the pseudo-inverse 
of r[2^. Inserting the parameterization of V and W in terms 
of Y defines a Hermitian x matrix, the eigenvalues of 
which produce variational estimates for the excitation ener- 
gies of the LL Hamiltonian. The lowest eigenvalues can then 
be obtained with a computational time scaling as O(D^) using 
a sparse eigensolver exploiting the tensor product structure of 
the expression above. 

We now discuss the results obtained for the LL model. 
Although the exact ground state does not break the U(l)- 
symmetry of the model, the cMPS approximation for the 
ground state does break this symmetry, just like the mean 
field solution (which corresponds to a cMPS with D = 1) 
obtained from the Gross-Pitaevskii equation. This follows 
from the fact that it is often energetically beneficial to break 
the symmetry in the presence of a constraint on the to- 
tal amount of entanglement. Hence, the order parameter 
('i'{Q,R)\-ip\'i'{Q,R)} / (and in fact slowly converges to 
for increasing D). For any 6 G [0, 2n), |*(Q, e'^i?)) is again 
a valid ground state so that we can also consider topological 
excitations interpolating between two different ground states 
characterized by a different order parameter. 

Fig. [T] shows the eigenvalues of the effective Hamiltonian 
Hp as a function of the momentum p for three different val- 
ues of 7 obtained using the trivial ansatz (Qi = Q2 and 
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FIG. 1. Dispersion relation for the Lieb-Liniger model with D = 
64 for 7 « 0.1,1.35,60.16 (up-down). Bulk excitations (dots) are 
trivial excitations. Blue diamonds (Type I) and red stars (Type II) 
are obtained by combining momenta and energy of two topological 
excitations (hole and particle) of Fig. [2] 



Ri = 7?2). All trivial eigenvalues are contained in a gap- 
less continuum. The elementary excitations in the LL model 
are the hole and the particle excitations which are studied in 
the Bethe ansatz by using anti-periodic boundary conditions 
[30l |. Fig. 2 shows that these excitations are perfectly cap- 
tured by using the ansatz for topological excitations with the 
choice Qi = Q2 ~ Q, Ri = —R2 = R- The dispersion relation 
is now centered around and the hole branch (red stars) is 
obtained as the lowest excitation energy between — tt and tt. 
The particle branch (blue diamonds) shows the eigenvalues of 
the eigenvectors that have the largest overlap with the state 
J^^ dx e'^'Vit (a:)e'" '^^("''^("' |*(Q, R)). 

Lieb determined the spectrum with fixed particle number 
in first quantization [2^ and isolated two types of excitations, 
which he labeled Type I and Type II excitations. Either can 
be used to construct the full spectrum of excitations with 
equal particle number. In our second quantized approach, we 
can understand the Type I excitations by combining the mo- 
mentum and energy of one hole at momentum —np plus one 
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FIG. 2. Elementary topological excitation branches for 7 = 60.16 
at D = 64 with 6 = tt. Blue diamonds are particle excitations and 
red stars hole excitations, (a) Dispersion relation of Lieb's Type 
I excitation obtained by combining energies from the topological 
ansatz (blue diamonds) and from the exact Bethe ansatz solution 
(black line), (b) Dispersion relation of Lieb's Type II excitation 
obtained from the topological ansatz (red stars) and from the exact 
Bethe ansatz solution (black line). The lowest eigenvalues of the 
trivial ansatz (green squares) are shown for comparison. 



particle with momentum p > np [blue diamonds in Fig. [T] and 
[2ja)]. The Type II excitations at momentum p are obtained 
by combining energy and momentum of one particle with mo- 
mentum np plus one hole with momentum —np < p < np [red 
stars in Fig. [T]and[2I^b)]. Using this construction, we can re- 
produce the Bethe ansatz result for the dispersion relation of 
Lieb's Type I and II excitations perfectly (5 digits of preci- 
sion), as indicated in Fig.[2Ia) and[5Jb). Due to this compos- 
ite structure, much less accurate energies would have resulted 
from using the trivial ansatz directly, as is also indicated in 
Fig. Ob). Since the repulsive Lieb-Liniger model does not 
have any bound states [30| , all trivial excitations have a simi- 
lar structure consisting of an even sum of topological hole and 
particle excitations. 

Without knowing the exact result, this information could 
be inferred from looking at the convergence of the variational 
energies as a function of increasing bond dimension. There are 
two competing effects within our variational strategy. Firstly, 
by modifying Q and R locally, we assume that the excita- 
tion is confined in a spatial region, the width of which is set 
by the bond dimension. This is the variational approxima- 
tion and it produces a positive variational error. Secondly, 
we subtract from the Hamiltonian a variational estimate of 
the ground state energy which is too large. This second effect 
results in a negative error and is dominant for truly elemen- 
tary excitations for which the assumed locality is valid and 
the first error is negligible. If an excited state cannot be 
approximated locally, e.g. an excitation that is composed of 
several unbound elementary excitations, then the first effect 
will likely be dominant as the ansatz confines these different 
excitations in a finite spatial region. Indeed, for increasing 
D, the particle and hole energies obtained from Fig. [2] are 
increasing (for 7 > 3) because a smaller variational estimate 
for the ground state energy is subtracted. In contrast, the 
energies of trivial excitations in Fig. [T] decrease for increasing 
bond dimension, in line with our expectations. 

In many ways, the topological hole excitation can be in- 
terpreted as a soliton. This observation was first discussed 
in Refs. [31I andls^ and recently received renewed attention in 
Ref. [s^ andls^ following the experimental realization [2^. In 
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FIG. 3. The change in particle number expectation value for the 
hole state at momentum p = as a function of 7, as predicted by 
the Bethe solution |34|| and with the cMPS excitation ansatz. 



the Tonks- Girardeau limit (large 7) [3^, the hole excitation 
maps exactly to the fermion-like hole excitation and produces 
a change in particle number of —1. For very small 7, the hole 
excitation is related to the classical dark soliton of the non- 
linear Schrodinger equation [3^ and produces a large change 
in particle number [3J]. The relative change in particle num- 
ber can easily be calculated using our ansatz. It is illustrated 
in Fig [3] at momentum p = 0, and the results perfectly coin- 
cide with the recent exact calculations [33 |. 

Next, we study a non-integrable model obtained by adding 
a pairing term to Hll that breaks the U(l) symmetry down 
to a residual Z2 symmetry 

H' = Hll + [ dx (uiP^'' +U'4iip) . 
J —00 

As long as v < the cMPS ground state approximation 
\'^(Q,R)) spontaneously breaks the Z2 symmetry for any 
nonzero pairing strength u £ Co, and a second ground 
state |vl/(Q,— i?)) is obtained. This also opens up a gap 
in the spectrum. The magnitude of the order parameter 
{'i!{Q, R)\^{x)\'i'{Q, R)) is determined by the competing ef- 
fects of the pairing term and the interaction term. As an ex- 
ample we consider an intermediate parameter range where we 
have found strong evidence that the lowest lying excitation in 
the trivial spectrum is a bound state. In the top plot of Fig. 01 
the lowest lying trivial excitations of H' for 7 « 26.4, v = —1 
and u — 1 with D — 22 are shown with red circles. Due to 
the symmetry breaking, we can again construct topological 
excitations (right insert of Fig. |4]|. The lowest lying eigen- 
values constitute an isolated branch that can be interpreted 
as the kink excitation that interpolates between the two de- 
generate ground states. The blue dots in the trivial spectrum 
of Fig. |4] are obtained by considering all possible pairs of two 
such topological excitations by adding their momentum and 
energy. Around momentum zero, the lowest lying trivial ex- 
citation produced by our topologically trivial ansatz lies well 
below the two-kink continuum starting at twice the kink mass. 

It is tempting to interpret this excitation as a bound state of 
two kinks, and this claim is further supported by considering 
the convergence behavior of (twice) the kink mass and the 
two lowest trivial excitation energies as a function of the bond 
dimension D in the bottom plot of Fig. |4] The kink is the truly 
elementary excitation; it has a negative error that vanishes 
quickly for increasing bond dimension. As in the LL case, 
unbound multi-particle states are not expected to converge 
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FIG. 4. Top: Dispersion relation of H' for 7 26.4, v = — 1 
and M = 1 . Blue dots are the pairwisc sum of all topological 
excitations. Red circles are trivial excitations and the right insert 
shows the topological spectrum. Bottom: Convergence of the two 
lowest trivial eigenvalues and the lowest two-kink excitation energy 
as a function of at p = 0. 



quickly, because the different constituents are confined within 
a spatial region, the size of which is determined by the bond 
dimension [37|l. For the energy of a two-particle state, we thus 
expect a positive finite size error. The energy of a bound state 
should also exhibit a large positive error as long as the spatial 
support of the excitation ansatz is below the binding length, 
but should then stagnate when it exceeds the binding length 
by increasing D. Fig. [4] confirms that the lowest lying energy 
of the trivial ansatz converges to a value well below twice the 
kink mass in a way that is characteristic for a bound state. 

In conclusion, we have introduced a variational ansatz for 
the study of elementary excitations for one-dimensional quan- 
tum field theories, based on the concept of continuous matrix 
product states. The ansatz produces very accurate results 
even for a critical model such as the interacting Bose gas. 
The local aspect of the ansatz made the solitonic feature of 
the lowest lying hole eigenstates of the Lieb-Liniger model 
very explicit, and the particle number expectation value in 
those solitons could be reproduced to great accuracy. We 
also applied the methods to a gapped non-integrable variant 
of LL, and obtained strong evidence for the existence of a 
bound state. 
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